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Abstract 

We introduce a lattice Boltzmann model for simulating an immiscible binary 
fluid mixture. Our collision rules are derived from a macroscopic thermo- 
dynamic description of the fluid in a way motivated by the Cahn-Hilliard 
approach to non-equilibrium dynamics and ensure that a thermodynamically 
consistent state is reached in equilibrium. The non-equilibrium dynamics is 
investigated numerically and found to agree with simple analytic predictions 
in both the one-phase and the two-phase region of the phase diagram. 

PACS numbers: 02.70Ns, 64.90. +b, 47.11. +j 
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Since the introduction of automaton based models capable of simulating hydrodynamic 
phenomena [I], lattice gas and lattice Boltzmann techniques have proved to be extremely 
useful numerical tools for investigating fluid dynamics. These approaches have been applied 
to a variety of different physical problems ranging from the simple flow of one component- 
fluids to multi-component fluid dynamics in random porous media [TJJ, a problem of great 
interest to the oil industry [[|. 

A feature common to all lattice Boltzmann schemes is to describe the fluid or fluid mix- 
ture by a set of distribution functions. These probabilities are assumed to obey a Boltzmann 
equation, discrete in both time and space, which is readily evolved numerically. The macro- 
scopic fluid variables, for example the density and velocity, are defined via moments of the 
distribution functions, and, provided the Boltzmann collision operator is chosen in a suitable 
way, can be shown to obey hydrodynamic transport equations ||]. 

For a single- component fluid with an ideal equation of state, the collision operator need 
only respect local conservation of mass and momentum. This is analogous to the Boltzmann 
approximation in continuum kinetic theory |5|]. However, non-ideal fluids and immiscible 
fluid mixtures are examples of interacting systems and to model the consequences of the 
interactions, most notably phase separation, effective collision processes which mimic the 
particle interactions must be defined. There have been several attempts to do this either 
by introducing effective interactions directly between the probability variables or by a 
phenomenological rewriting of the collision rules to favour interface formation 0. Such 
approaches have proved successful in a variety of applications ||, but all have the disad- 
vantage that the resulting macroscopic behaviour cannot, in general, be related directly to 
a thermodynamic description of the fluid 0. 

Therefore in this letter we introduce a lattice Boltzmann scheme for simulating immisci- 
ble binary liquid mixtures which retains the numerical efficiency of earlier lattice Boltzmann 
models of multi-phase systems, while ensuring the equilibrium macroscopic behaviour is 



consistent with a thermodynamic description of the fluid mixture flIU|. Our approach was 



motivated by considering the lattice Boltzmann technique as a mesoscopic rather than mi- 



croscopic description of a fluid which allowed us to input the conventional coarse-grained 
description of non-equilibrium dynamics |TT |. 



The basic physical variables we shall choose to work with are (1) the total fluid density 
p = pi + P2, (2) the mean fluid velocity u and (3) the density difference between the two 
components Ap = pi — P2, where p± and pi are the individual component densities. We 
introduce two distribution functions each of which evolves according to a single relaxation 



time Boltzmann equation [12 



fi(x + e t At,t + At) - n{x,t) = --{h - f^), (1) 



A i (x + e i At,t + At) - Ai(x,t) = — — (Aj - Af), (2) 

ta 

where t p and ta are independent relaxation parameters and is a lattice vector. The 
distribution functions are related to the physical variables by 

P = J2fi> PU=J2fi^ Ap = J2^i. (3) 

i i i 

These three quantities are locally conserved in any collision process and this requirement im- 
poses three constraints on the equilibrium distribution functions in eqns(l) and (2), namely, 

E/r = P, E/fe^pM, £Af = Ap. (4) 

i i i 

The higher moments of fl q and A eq are defined so that the resulting continuum equations 
take the form pertinent to a binary liquid mixture [O]. We require that Ap obeys the Cahn- 
Hilliard equation supplemented by an advection term, while the fluid as a whole obeys the 
Navier-Stokes equations with a non-ideal pressure tensor. To this end we define 

£ Afe ioi = Apu a , e ia e i/3 = TApS af3 + Apu a u p , (5) 

i i 

for the second and third moment of A^ 9 , and 

fi qe ia e iP = Pa/3 + PU a Up (6) 



for the third moment of fl q . In the above equations Ap is the chemical potential difference 
between the two components, T is the mobility and P Qj g is the pressure tensor. 

With these definitions of f^ Q and A^ the continuum equations which follow from ex- 
panding eqns(l) and (2) to 0(At 2 ), are 

d t A P + d a A P u a = rev 2 A/2 - ed^dpPnp, (7) 

dtp + d a pu a = 0, (8) 



d t pu a + dppuaup = -dpP al3 + uV 2 pu a + 4ud a (XV.pu), (9) 

where 

= A(c* (r A - i) . v = 1*»\{t,-\), A = g - g . (10) 

Note that the viscosity is controlled by r p while the introduction of the coefficient T allows 
the diffusion constant to be varied independently of ta- 

If p is constant eqns(7-9) are Galilean invariant. However, there will be correction terms 
of 0(At 3 ), and higher, which break this invariance. By expanding eqn(2) to third order in 
At we find a ta dependent pre-factor P^(ta) = Ta — r A + |- Thus, by choosing ta to be a 
zero of P^(ta) we can eliminate all third order corrections to the diffusion equation while 
retaining control over the diffusion constant through T. Numerically, we find this to be an 
important requirement in order to preserve Galilean invariance as closely as possible. 

The thermodynamic aspects of the model enter through Ap and P a p. Following the 
Cahn-Hilliard description of non-equilibrium dynamics, we calculate both of these functions 
from the equilibrium free energy of the fluid mixture. We choose the simplest model of a 
binary liquid: two ideal gases with a repulsive interaction energy. In terms of our model 
variables, the free energy functional of the mixture takes the form 

* = / d 2 r{^(Ap,p,T) + |(V P ) 2 + ^(VAp) 2 }, (11) 
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in which ip(Ap, p, T) is the bulk free energy density at a temperature T and the second two 



terms give the free energy contribution from density gradients ||13|| . For two interacting ideal 

gases 



\p(l - ^) - Tp + |(p + Ap) log l e+^L) + T - (p - Ap) log (Zl££) 



where A measures the strength of the interaction. For T < T c = |A the bulk system phase 
separates into one of two phases, symmetric in Ap. From this free energy the chemical 



potential difference and pressure tensor follow in the usual way [13 



Ap(Ap,p,T) = -A^ + Tlog ~ «V 2 (Ap), (13) 



with 



p(r) = P T-\ (pV 2 p + ApV 2 Ap) - | (|Vp| 2 + | VAp| 2 ) . (15) 

The distribution functions f* 9 and Af 3 must be chosen to satisfiy eqns(4-6), which in- 
corporate the thermodynamic description of the fluid. As only the first three moments of 
fl q and A eq are specified, it is sufficient to expand these functions to second order in u. 
The coefficients in the expansion are allowed to be not only functions of p and Ap but 
also functions of their derivatives. This ensures that Ap and P a p have the correct form in 
inhomogenous regions of the fluid. 

We now turn to the numerical implementation of the model in order to test the theory 
presented above and to establish a practical limit on the scheme. Most of the simulations 
were performed on a triangular lattice with 128 x 128 sites. The strength of the interaction 
A was set to 1.1 giving a critical temperature T c = 0.55. We chose k = 0.1, c = 1 and worked 
with units in which At = 1. The behaviour of the system is controlled by the temperature 
parameter T and we have investigated the equilibrium and dynamical properties in three 
distinct temperature regimes: T > T C ,T = T c and T < T c . 



Above the critical temperature the equilibrium configuration is a homogeneous mixture 
of the two components, characterised by p = constant and Ap = 0. For small variations 
in Ap, eqn(7) can be linearised about Ap = resulting in a convective-diffusion equation 
with diffusion coefficient D = 26T(T — T c ). To test this prediction we have measured D as 
a function of T by monitoring the decay of a sinusoidal perturbation in Ap. Figure 1 shows 
the measured values of D as a function of T for different values of T > T c . Agreement with 
the predicted dependence is excellent. We have also investigated the Galilean invariance 
inherent in the model by imposing a uniform velocity onto the initial perturbation. We 
find the measured values of D agree with the static case for a wide range of velocities: D 
varies by less than |% for < u < 0.5. Note that our approach in the miscible regime is 
an extension of the work by Flekkoy ]14| to include additional terms in A^ 9 which restore 



Galilean invariance. 

Strictly at the critical temperature, the diffusion coefficient discussed above is zero. We 
would thus expect a difference in the dynamics at T c even though the equilibrium state is 
uniform. This difference can be investigated by monitoring the decay of the non-equilibrium 
surface tension, defined for a single interface by a = J (^r) dx, where the integral is taken 
perpendicular to the gradient in Ap. Two distinct power laws are observed for the decay 
of a, depending on whether the final temperature Tf is above or equal to T c : a ~ t~^ for 
Tj > T c and o ~ t~^ for Tf = T c . This is in agreement with the predictions of Ma et. al. 
|T5| and indicates that the system is in the Model B universality class, in a regime where 
hydrodynamics is unimportant. Earlier lattice Boltzmann models of binary mixtures were 
unable to reproduce such details. 

Finally, below T c the fluid phase separates into two distinct phases, symmetric about 
Ap = 0. Figure 2 shows typical interfacial profiles between two phase characterised by 
Ap > and Ap < 0. The width of the interface can be tuned by varying T (or the 
interfacial energy k) and chosen so that lattice effects are unimportant. We also show the 
variation of the total density p in the interfacial region; typically, this variation is less then 
2%. Far from the interfacial region, the uniform values of Ap define the coexisting phases. 



The coexistence curve as a function of T is shown in figure 3. The observed values of Ap 
are seen to be in agreement with the theoretical prediction derived by minimising the bulk 
free energy, eqn(12). 

We have verified that deviations from Galilean invariance are very small even in the 
immiscible regime: there is little change to either the coexisting values of Ap or the structure 
of the interface if a uniform velocity is imposed on the system. The velocity dependence of the 
coexistence curve is also show in figure 3. A well know problem common to all multi-phase 
simulation techniques is the generation of spurious velocity currents in the interfacial region. 
We have measured this velocity field within our model and find its magnitude appears to be 
substantially lower than the values that have been observed using other schemes J?||10| . For 
example, for the interfaces shown in figure 2, \u\ is O(10~ 10 ). Lastly, we have checked that 
Laplace's Law holds for spherical domains, a feature which is ensured in the continuum by 
our choice of pressure tensor. A detailed study of the collective dynamics during spinodal 
decomposition will be presented elsewhere |LJ|. 

To summarise: we have devised and tested a lattice Boltzmann scheme for simulating 
a binary liquid mixture which leads to a thermodynamically consistent equilibrium state. 
Three distinct dynamical regimes are observed; above, at and below the critical temperature, 
and the agreement between our simulation results and theoretical predictions is excellent 
for a wide range of system parameters. Furthermore, the effects of correction terms which 
break Galilean invariance are seen to be small, even in the presence of sharp interfaces. 

We are indebted to William Osborn and Jayanth Banavar for numerous valuable discus- 
sions. EO acknowledges support from the European community and MRS and JMY from 
the E. P. S. R. C. 
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FIGURE CAPTIONS 



Figure 1 The diffusion coefficient D as a function of the parameter T for three different 
values of the temperature. The points represent results from simulations: open circles 
T=0.8, filled square T=0.7 and filled circles T=0.6. The solid lines show the analytic 
value D = 2T6 (T-T c ). 

Figure 2 Profiles normal to a flat interface of the binary mixture for three different values 
of the temperature: T=0.498 (solid), T=0.511 (dotted) and T=0.526 (dashed). The 
parameter k is fixed to the value 0.1. Note that the variation of p in the interfacial 
region is small; the 'spikes' around x = are an artifact of lattice discretisation. 

Figure 3 Coexistence curve as a function of the temperature T for A = 1.1. The different 
points denote simulations at different uniform fluid velocity: filled squares u=0.0, open 
triangles u=0.1 and open hexagons u=0.2. The solid line is the analytic curve obtained 
by minimising the bulk free energy. 



10 



0.05 



0.04 



0.03 



0.02 



0.01 







